#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <math.h>
#include <string.h>
#include "mkl_solvers_ee.h"
#include <fstream>
#include "mkl.h" // for mkl_malloc and mkl_free
#define max(a, b) (a) < (b) ? (b): (a)
#define min(a, b) (a) < (b) ? (a): (b)
#define ALIGNMENT 128  //4096
int main()
{
    char UPLO = 'F';
    int nrows=2, nnz=2;
    

    MKL_INT *rowPtr=NULL,*colInd=NULL;
    double *values=NULL;
    rowPtr = (MKL_INT*)mkl_malloc((nrows+1)*sizeof(MKL_INT), ALIGNMENT);
    colInd = (MKL_INT*)mkl_malloc(nnz*sizeof(MKL_INT), ALIGNMENT);
    values =  (double*)mkl_malloc(nnz*sizeof(double), ALIGNMENT);
    // 非零元素的值
    // std::vector<double> values = {3, 2};
    // // 非零元素的列索引
    // std::vector<MKL_INT> colInd = {0, 1};
    // // 每行非零元素的起始位置
    // std::vector<MKL_INT> rowPtr = {0, 1, 2};


    for(unsigned int k=0;k<nnz;k++) values[k]=-1.0;// set values to -1

    /* Declaration of FEAST variables */
    MKL_INT fpm[128];    /* Array to pass parameters to Intel MKL Extended Eigensolvers */
    double  Emin, Emax;  /* Lower/upper bound of search interval [Emin,Emax] */
    double  epsout;      /* Relative error on the trace */
    MKL_INT loop;        /* Number of refinement loop */
    MKL_INT L = 2;       /* expected eigenvalues to be found */
    MKL_INT M0;          /* Initial guess for subspace dimension to be used, should be larger than L */
    MKL_INT M;           /* Total number of eigenvalues found in the interval */
    double  *E=NULL;     /* Eigenvalues */
    double  *X=NULL;     /* Eigenvectors */
    double  *res=NULL;   /* Residual */

    /* Declaration of local variables */
    MKL_INT      info;          /* Errors */
    double       one = 1.0;     /* alpha parameter for GEMM */
    double       zero = 0.0;    /* beta  parameter for GEMM */
    MKL_INT      i, j;
    double       trace, smax, eigabs;

    printf("\n FEAST DFEAST_SCSREV AND DFEAST_SCSRGV EXAMPLE\n"); fflush(0);

    /* Search interval [Emin,Emax] */
    Emin = 0; //smallest eig is ~ -12.75531,  largest is ~12.5531,  they are all very close together in this interval.  You can't give whole interval since that is 93k eigs... too much for dense subspace 
    Emax = 4;
    printf("Search interval [ %.15e, %.15e  ]  \n", Emin, Emax); fflush(0);

    M0   = min(L+L/4, nrows);  // set M0 a little larger than expected number of eigs to be found... sometimes it needs to be even larger...
    M    = L;                  // expected number of eigs to be found
    loop = 3;
    info = 0;
    epsout = 0.0;

    // set arrays to proper size  should be room for subspace size M0 eigenvalues and eigenvectors
    E   = (double*)mkl_malloc(M0*sizeof(double),ALIGNMENT);
    X   = (double*)mkl_malloc(M0*nrows*sizeof(double),ALIGNMENT);
    res = (double*)mkl_malloc(M0*sizeof(double),ALIGNMENT);

    /* Initialize matrix X */
    for (i=0; i<L*nrows; i++){ // nrows*nrows
        X[i] = zero;
    }

     /* Step 1. Call  FEASTINIT to define the default values for the input FEAST parameters */
    feastinit(
        fpm /* OUT: Array is used to pass parameters to Intel MKL Extended Eigensolvers */
        );

    fpm[0] =  1; /* Extended Eigensolver routines print runtime status to the screen. */

    /* Step 2. Solve the standard Ax = ex eigenvalue problem. */
    printf("Testing dfeast_scsrev routine:\n"); fflush(0);

    dfeast_scsrev(
            &UPLO,   /* IN: UPLO = 'F', stores the full matrix */
            &nrows,      /* IN: Size of the problem */
            values,     /* IN: CSR matrix A, values of non-zero elements */
            rowPtr,    /* IN: CSR matrix A, index of the first non-zero element in row */
            colInd,    /* IN: CSR matrix A, columns indices for each non-zero element */
            fpm,     /* IN/OUT: Array is used to pass parameters to Intel MKL Extended Eigensolvers */
            &epsout, /* OUT: Relative error of on the trace */
            &loop,   /* OUT: Contains the number of refinement loop executed */
            &Emin,   /* IN: Lower bound of search interval */
            &Emax,   /* IN: Upper bound of search interval */
            &M0,     /* IN: The initial guess for subspace dimension to be used. */
            E,       /* OUT: The first M entries of Eigenvalues */
            X,       /* IN/OUT: The first M entries of Eigenvectors */
            &M,      /* OUT: The total number of eigenvalues found in the interval */
            res,     /* OUT: The first M components contain the relative residual vector */
            &info    /* OUT: Error code */
            );

    printf("FEAST OUTPUT INFO %d \n",info); fflush(0);
    if ( info != 0 )
    {
        printf("Routine dfeast_scsrev returns code of ERROR: %i", (int)info);
        return 1;
    }

    printf("Number of eigenvalues found %d \n", M); fflush(0);
    for (i=0; i<M; i++){
        printf("%.15e \n", E[i]);
    }
    
    
release_memory:
    mkl_free(rowPtr); mkl_free(colInd); mkl_free(values);
    mkl_free(E); mkl_free(X); mkl_free(res);

    return 0;

}